N2O AF calcs v2¶

12/06/24

  • Rerun for trange=[17,23]
  • Use updated HV plotting routines.

22/03/24

Alignment frame calcs preliminary:

  • Load and plot alignment data (ADMs) for all temperatures.
  • Compute AF results for given temp (with t downsampling), note default case assumes linear pol, $z$-axis alignment.
  • Plot XS and betas, interactive heatmaps and line plots.

Note all plots are interactive, mouse-over for data points, and use toolbars for zoom/selection, and widgets for browsing data types (where applicable).

For general methods: https://epsproc.readthedocs.io/en/dev/demos/ePSproc_class_demo_161020.html

In [1]:
# Quick hack to override default HTML template
# NOT required in some JLab versions.
# https://stackoverflow.com/a/63777508
# https://stackoverflow.com/questions/21971449/how-do-i-increase-the-cell-width-of-the-jupyter-ipython-notebook-in-my-browser
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))

/tmp/ipykernel_18245/28059362.py:5: DeprecationWarning: Importing display from IPython.core.display is deprecated since IPython 7.14, please import from IPython display
  from IPython.core.display import display, HTML

Load alignment data¶

Now with ADM class.

In [2]:
from epsproc.classes.alignment import ADM
from pathlib import Path

# dataPath = Path('~/ePS/N2O/alignment/Wigner D polyatomic normalized')
dataPath = Path('~/jake-home/ePS/N2O/alignment/Wigner D polyatomic normalized')
ADMtype = '.txt'

ADMin = ADM(fileBase=dataPath.expanduser(), ext = ADMtype)

ADMin.loadData(keyType='setADMsT', normType='wignerDpoly')

OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.

* sparse not found, sparse matrix forms not available. 
* natsort not found, some sorting functions not available. 

* Setting plotter defaults with epsproc.basicPlotters.setPlotters(). Run directly to modify, or change options in local env.

* Set Holoviews with bokeh.
* pyevtk not found, VTK export not available. 
Scanning /home/jovyan/jake-home/ePS/N2O/alignment/Wigner D polyatomic normalized for ePS jobs.
Found ePS output files in subdirs: [PosixPath('/home/jovyan/jake-home/ePS/N2O/alignment/Wigner D polyatomic normalized/A200'), PosixPath('/home/jovyan/jake-home/ePS/N2O/alignment/Wigner D polyatomic normalized/A400'), PosixPath('/home/jovyan/jake-home/ePS/N2O/alignment/Wigner D polyatomic normalized/A600'), PosixPath('/home/jovyan/jake-home/ePS/N2O/alignment/Wigner D polyatomic normalized/A800')]

*** Scanning dir
/home/jovyan/jake-home/ePS/N2O/alignment/Wigner D polyatomic normalized/A200
Found 6 .txt file(s)


*** Scanning dir
/home/jovyan/jake-home/ePS/N2O/alignment/Wigner D polyatomic normalized/A400
Found 6 .txt file(s)


*** Scanning dir
/home/jovyan/jake-home/ePS/N2O/alignment/Wigner D polyatomic normalized/A600
Found 6 .txt file(s)


*** Scanning dir
/home/jovyan/jake-home/ePS/N2O/alignment/Wigner D polyatomic normalized/A800
Found 6 .txt file(s)

Set self.norm from self.norms['wignerDpoly'].
Found t-index /home/jovyan/jake-home/ePS/N2O/alignment/time.txt
Files from self.fileDict read OK, outputs: self.dataDict and self.data.
In [3]:
# import epsproc as ep

# ep.setADMs
In [4]:
# ADMin.plot(width=1000)
# ADMin.plot(xlim=(30,50), width=800)  # Auto-stack plus pass args & zoom in on specific region (note slider will reset region with interactive zoom)

Load matrix elements¶

In [5]:
# Plotters
# from epsproc.plot import hvPlotters

# Multijob class dev code
from epsproc.classes.multiJob import ePSmultiJob

import warnings
# warnings.filterwarnings('once')   # Skip repeated numpy deprecation warnings in current build (xr15 env)
warnings.filterwarnings('ignore')   # Skip repeated numpy deprecation warnings in current build (xr15 env)
In [6]:
# # Scan for subdirs, based on existing routine in getFiles()

# fileBase = Path('/home/paul/ePS/OCS/OCS_survey')  # Data dir on Stimpy

fileBase = Path('~/fock-mount/globalhome/eps/N2O/N2O_valence').expanduser() # Data dir on Jake
In [7]:
# TODO: fix orb label here, currently relies on (different) fixed format

data = ePSmultiJob(fileBase, verbose = 0)

data.scanFiles()
# data.jobsSummary()

*** Warning: Missing records, expected 32, found 24.
*** Warning: Found 8 blank sets of matrix elements, symmetries ['A2']
*** Warning: Missing records, expected 32, found 24.
*** Warning: Found 8 blank sets of matrix elements, symmetries ['A2']
*** Warning: Missing records, expected 32, found 24.
*** Warning: Found 8 blank sets of matrix elements, symmetries ['A2']
*** Warning: Missing records, expected 32, found 24.
*** Warning: Found 8 blank sets of matrix elements, symmetries ['A2']
*** Warning: Missing records, expected 32, found 24.
*** Warning: Found 8 blank sets of matrix elements, symmetries ['A2']
*** Warning: Missing records, expected 32, found 24.
*** Warning: Found 8 blank sets of matrix elements, symmetries ['A2']
*** Warning: Missing records, expected 32, found 24.
*** Warning: Found 8 blank sets of matrix elements, symmetries ['A2']
*** Warning: Missing records, expected 32, found 24.
*** Warning: Found 8 blank sets of matrix elements, symmetries ['A2']

Compute AF¶

Note default case assumes linear pol, $z$-axis alignment.

In [8]:
# Set ADMs to use for calc - set for specific T and downsample

# Subselection, note temperature is set by dataKey here
temp = '10K'
trange = [17,23]
# ADMin.subsetADMs(dataKey = temp, trange=[35,45],tStep = 5)
ADMin.subsetADMs(dataKey = temp, trange=trange,tStep = 2)

# Plot subselection
ADMin.plot(keys='ADM')

# Set to main data structure for calcs
data.data['ADM'] = ADMin.data['ADM']

Setting subset data from `self.data['10K']['ADM']`
Selecting 60 points
Set subset data to `self.data['ADM']['ADM']`
In [9]:
# # Basic routine (no grouping)
# import epsproc as ep
# thres = 1e-5
# # data.AFBLM(AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=1e-2) # OK for 41 t-points (x2 downsample), but missing some orbs and t-points in output!
# # data.AFBLM(AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=thres) # OK for 41 t-points (x2 downsample), still some missing regions for thres=1e-2
#                                                                                        # mostly contiguous for thres = 1e-3,1e-4
# data.AFBLM(AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=thres)   # outputs fairly smooth for thres=1e-5

# # Save full class with Pickle
# # TODO: GENERAL DEBUG FOR ACCIDENTAL FUNCTION PASSING HERE!!!!
# import pickle
# outfile = f'N2O_AFBLM_{temp}_t{trange[0]}-{trange[1]}_{thres}_120624.pickle'
# with open(outfile, 'wb') as handle:
# #     data.data['ADM']['plots']['ADM']['pData'].attrs['pTypeDetails'] = []  # UGH, this breaks pickle unless removed.
#     pickle.dump(data.data, handle, protocol=pickle.HIGHEST_PROTOCOL)
#     print(f'Saved data.data with pickle, file: {outfile}')
    
#     # Also save to HDF5 (AFBLM only)
#     for orb in data.data.keys():
#         # Skip non-orb keys
#         if orb.startswith('orb'):
#             ep.IO.writeXarray(data.data[orb]['AFBLM'], fileName=f'N2O_AFBLM_{orb}_{temp}_t{trange[0]}-{trange[1]}_{thres}_120624.h5', engine='hdf5')
In [10]:
# Load data from previous runs...
thres = 1e-5
# infile = f'N2O_AFBLM_{temp}_{thres}_220323.pickle'
# infile = 'N2O_AFBLM_10K_0.0001_220323.pickle'
# infile = 'N2O_AFBLM_10K_1e-05_220323.pickle'
infile = f'N2O_AFBLM_{temp}_t{trange[0]}-{trange[1]}_{thres}_120624.pickle'

import pickle
with open(infile, 'rb') as handle:
    dataIn = pickle.load(handle)
    
data.data = dataIn

Plots¶

Remarks:

  • Some unphysical spikes in values <5eV, these are regions with small cross-sections, and are currently filtered out in these plots.
  • General patterns seem to be simlar to R-matrix results, as expected, although absolute Blm values seem larger...? Might be a normalisation difference here, TBC.

Cross-sections & AF-$\beta_{L,M}$ heatmaps¶

NOTE: interactive plots. Roll over plots for values and tool options. For heatmaps use the "box select" tool on the histogram at the side to change the colour mapping scale.

In [11]:
# v2 with pkg code...
# from epsproc.plot import hvPlotters
# hvPlotters.setPlotDefaults(fSize = [700,300], imgSize = 600)
# data.BLMplot(backend='hv', xDim = 't', hvType='heatmap', addHist = True, addADMs = True, XS=True, 
#              clim=(-1,2), cmap='vlag', width=1000, frame_width=860, filterXS=1e-3)
In [12]:
data.data['orb8_S']['AFBLM'].t[slice(20,40,5)]    #[5,10]]    #slice(5,10,2)]
Out[12]:
<xarray.DataArray 't' (t: 4)>
array([19.022921, 19.52086 , 20.018799, 20.516738])
Coordinates:
  * t        (t) float64 19.02 19.52 20.02 20.52
    Type     <U1 'L'
    it       int64 1
Attributes:
    units:    ps
xarray.DataArray
't'
  • t: 4
  • 19.02 19.52 20.02 20.52
    array([19.022921, 19.52086 , 20.018799, 20.516738])
    • t
      (t)
      float64
      19.02 19.52 20.02 20.52
      units :
      ps
      array([19.022921, 19.52086 , 20.018799, 20.516738])
    • Type
      ()
      <U1
      'L'
      array('L', dtype='<U1')
    • it
      ()
      int64
      1
      array(1)
  • units :
    ps

AFPADs¶

Selected by orb, Eke and t.

In [13]:
# data.verbose['main']=True
# data.padPlot(keys='orb8_S', dataType='AFBLM',Etype='t', squeeze=True, selDims={'Eke':[5.1],'t':data.data['orb8_S']['AFBLM'].t[slice(20,50,5)]}, 
#              width = 1400, height = 800, backend='pl')
In [14]:
# AFPADs at selected (E,t)
# data.verbose['main']=True
data.padPlot(keys='orb8_S', dataType='AFBLM',Etype='t', squeeze=True, selDims={'Eke':[5.1],'t':data.data['orb8_S']['AFBLM'].t[slice(20,50,2)]}, 
             width = 1400, height = 1200, backend='pl')

Using default sph betas.
Summing over dims: set()
Plotting from self.data[orb8_S][AFBLM], facetDims=['t', None], pType=a with backend=pl.
*** WARNING: plot dataset has min value < 0, min = (-0.10478003837370972-8.928897465904397e-18j). This may be unphysical and/or result in plotting issues.
In [15]:
# Phi-summed PADs vs. t (single Eke)
data.padPlot(dataType='AFBLM',Etype='t', squeeze=False, selDims={'Eke':[5.1]},   # selDims={'t':data.data['orb8_S']['AFBLM'].t[slice(20,50,2)]}, 
             width = 1400, height = 500, reducePhi='sum', pStyle='grid', backend='hv')  #, hvClims=0.5,)  # facetDims=['Eke',None])  #, hvClims=None, hvCHist=False)   #facetDims=['Labels', 'Eke'])  #, backend='pl', reducePhi='sum')

Using default sph betas.
Summing over dims: set()
Found additional dims {'Eke'}, summing to reduce for plot. Pass selDims to avoid.
Plotting from self.data[orb6_S][AFBLM], facetDims=['t', 'Labels'], pType=a with backend=hv.
Grid plot: orb6_S, dataType: AFBLM, plotType: a
Using default sph betas.
Summing over dims: set()
Found additional dims {'Eke'}, summing to reduce for plot. Pass selDims to avoid.
Plotting from self.data[orb7_P][AFBLM], facetDims=['t', 'Labels'], pType=a with backend=hv.
Grid plot: orb7_P, dataType: AFBLM, plotType: a
Using default sph betas.

WARNING:param.HeatMapPlot02201: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.
WARNING:param.HeatMapPlot02201: aspect value was ignored because absolute width and height values were provided. Either supply explicit frame_width and frame_height to achieve desired aspect OR supply a combination of width or height and an aspect value.

Summing over dims: set()
Found additional dims {'Eke'}, summing to reduce for plot. Pass selDims to avoid.
Plotting from self.data[orb8_S][AFBLM], facetDims=['t', 'Labels'], pType=a with backend=hv.
Grid plot: orb8_S, dataType: AFBLM, plotType: a
Using default sph betas.
Summing over dims: set()
Found additional dims {'Eke'}, summing to reduce for plot. Pass selDims to avoid.
Plotting from self.data[orb9_P][AFBLM], facetDims=['t', 'Labels'], pType=a with backend=hv.
Grid plot: orb9_P, dataType: AFBLM, plotType: a

WARNING:param.HeatMapPlot02208: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.
WARNING:param.HeatMapPlot02208: aspect value was ignored because absolute width and height values were provided. Either supply explicit frame_width and frame_height to achieve desired aspect OR supply a combination of width or height and an aspect value.
WARNING:param.HeatMapPlot02213: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.
WARNING:param.HeatMapPlot02213: aspect value was ignored because absolute width and height values were provided. Either supply explicit frame_width and frame_height to achieve desired aspect OR supply a combination of width or height and an aspect value.
WARNING:param.HeatMapPlot02218: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.
WARNING:param.HeatMapPlot02218: aspect value was ignored because absolute width and height values were provided. Either supply explicit frame_width and frame_height to achieve desired aspect OR supply a combination of width or height and an aspect value.

AF-$\beta_{L,M}$ line-outs¶

Results as above, line-outs per (orb,Eke) with interactive browser widget.

In [16]:
# Line-outs
# Note this is quite slow to render
data.BLMplot(backend='hv', hvType='line', xDim = 't', ylim=(-1,2), filterXS=1e-3, addADMs = False, height=800, width=1000)

BLMplot set data and plots to self.plots['BLMplot']
Out[16]:
True

Versions¶

Code¶

In [19]:
# print(data.jobInfo['ePolyScat'][0])
# print('Run: ' + jobInfo['Starting'][0].split('at')[1])
In [20]:
import scooby
scooby.Report(additional=['epsproc', 'xarray', 'jupyter', 'holoviews', 'pandas'])
Out[20]:
Wed Jun 12 15:07:55 2024 EDT
OS Linux CPU(s) 64 Machine x86_64 Architecture 64bit
RAM 62.8 GiB Environment Jupyter File system btrfs
Python 3.10.11 | packaged by conda-forge | (main, May 10 2023, 18:58:44) [GCC 11.3.0]
epsproc 1.3.2-dev xarray 2022.3.0 jupyter Version unknown holoviews 1.16.2
pandas 1.5.3 numpy 1.23.5 scipy 1.10.1 IPython 8.13.2
matplotlib 3.5.3 scooby 0.7.2
In [21]:
# Check current Git commit for local ePSproc version
!git -C {Path(ep.__file__).parent} branch
!git -C {Path(ep.__file__).parent} log --format="%H" -n 1

/bin/bash: -c: line 1: syntax error near unexpected token `('
/bin/bash: -c: line 1: `git -C {Path(ep.__file__).parent} branch'
/bin/bash: -c: line 1: syntax error near unexpected token `('
/bin/bash: -c: line 1: `git -C {Path(ep.__file__).parent} log --format="%H" -n 1'
In [22]:
# Check current remote commits
!git ls-remote --heads https://github.com/phockett/ePSproc

1f115ae15db4357a87a7d6c9266b0285a94ecf52	refs/heads/3d-AFPAD-dev
18901ac6056324d6a2674396939fdcc373d35395	refs/heads/dependabot/pip/notes/envs/envs-versioned/idna-3.7
9d9eb4e1a20d4ff77d07af2ac782ba18c26e5dcc	refs/heads/dependabot/pip/notes/envs/envs-versioned/jinja2-3.1.4
226759a58ebe96bf1a6df60ccd5efb0c449af3a7	refs/heads/dependabot/pip/notes/envs/envs-versioned/pillow-10.3.0
de0e271701949602ce7f170a9665e37b27d2401c	refs/heads/dependabot/pip/notes/envs/envs-versioned/requests-2.32.0
4ff0ed84a1df2a3de4258ba7e49db1e47e6ac596	refs/heads/dependabot/pip/notes/envs/envs-versioned/scipy-1.11.1
120416a432b8da85b68912d0c424f47d77392434	refs/heads/dependabot/pip/notes/envs/envs-versioned/tornado-6.4.1
b50b0b946a1a7cb5b22c95d1e4cee120b22e6874	refs/heads/dependabot/pip/notes/envs/envs-versioned/tqdm-4.66.3
fd9c621f0c91e05ffcb097f59e9d8d8b43c5fc23	refs/heads/dependabot/pip/scipy-1.11.1
7e4270370d66df44c334675ac487c87d702408da	refs/heads/dev
1c0b8fd409648f07c85f4f20628b5ea7627e0c4e	refs/heads/master
69cd89ce5bc0ad6d465a4bd8df6fba15d3fd1aee	refs/heads/numba-tests
ea30878c842f09d525fbf39fa269fa2302a13b57	refs/heads/revert-9-master
baf0be0c962e8ab3c3df57c8f70f0e939f99cbd7	refs/heads/testDev

Notes¶

22/03/24

v1, quick run based on previous OCS code plus a few recent function updates.

TODO:

  • Testing with alternative polarization states.
  • More careful tests/checks in low-XS regions (currently just filter out, but see BLM spikes in a few cases without XS filter applied).
  • More careful comparison with R-matrix results.
  • AF-PAD plots, grids should work here...

Scratch¶

In [23]:
break

  Cell In[23], line 1
    break
    ^
SyntaxError: 'break' outside loop
In [ ]:
# Compute AFBLMs for aligned case
keys = data._keysCheck(None)
keys = list(set(keys) - set(['ADM']))


#*************** Throw everything in method - breaks for large datasets (RAM issues)
# data.AFBLM(keys = keys, AKQS = ADMs, selDims = {'Type':'L'})  # NOTE THIS USES independently set ADMs!

# data.AFBLM(keys = keys, AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=1e-3) # OK for 61 t-points (x2 downsample)
# data.AFBLM(keys = keys, AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=1e-3, thresDims = ['Eke','t'])   # NEED TO ADD DIM CHECKING TO thresDims
# data.AFBLM(keys = keys, AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=1e-4)   # 1e-4 working OK for 30 t-points (x4 downsample)
# data.AFBLM(keys = keys, AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=1e-5)   # 1e-5 working OK for 16 t-points (x8 downsample)
# data.AFBLM(keys = keys, AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=1e-6)   # 1e-6 working OK for 16 t-points (x8 downsample)
# data.AFBLM(keys = keys, AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=None)   # FAILS for 16 t-points (x8 downsample) at 60Gb RAM on Jake (Swap file size issue????)
                                                                                                   # OK for 8 t-points (x16 downsample) at ~55Gb RAM peak on Jake.
# thres=1e-4)  BREAKS KERNEL #, thresDims=['Eke','t'])  # TODO: currently thersDims must be in matE only for AFBLM routine.

# 21/07/22 version - above only gives L=2 max, missing L=4,6. But bumping threshold to 1e-4 kills calculation...
# data.AFBLM(keys = keys,  selDims = {'Type':'L'})  # Default case... sets/uses ADMs per key
# UPDATE: testing with scale factor for ADMs, and forcing to real - might be better?


#*********** Using grouping to avoid memory issues 
#        QUITE slow with current code, but more robust
#        TODO: make this better with basis set recycle (cf. fitting code)
#
import pickle

def calcAFBLMGroup(data,ADMs,keys,AFBLMdict):
    """Wrap for XR groupby"""
    
    # Class routine - returns to main datastruture
    data.AFBLM(keys = keys, AKQS = ADMs, selDims = {'Type':'L'}, thres=None)  #.map(standardize)
    
    # Copy to new structure & restack
#     keys = data._keysCheck(keys)

#     AFBLMdict = {}
    # Restack data to lists for merge
    for key in keys:
        AFBLMdict[key]['AFBLMlist'].append(data.data[key]['AFBLM'].copy())
        
    return AFBLMdict
        

AFBLMdict = {k:{'AFBLMlist':[]} for k in keys}
    
# tBins = np.arange(38,39.1,0.5)   # Test range
tBins = np.arange(trange[0], trange[1] + 0.1, 0.5)
# trange

import time

for item in data.data['ADM']['ADM'].groupby_bins("t", tBins):
    start = time.time()
    print(f'Calculating for group: {item[0]}, {item[1].t.size} t-points.')
    AFBLMdict = calcAFBLMGroup(data, item[1], keys, AFBLMdict)
    end = time.time()
    print(f'Delta: {end - start}s')

# Concat XR lists
for key in keys:
    AFBLMdict[key]['full'] = xr.concat(AFBLMdict[key]['AFBLMlist'],"t")
    AFBLMdict[key]['full'].name = f'{key}_AFBLM'
    
    # Push back to main class
    data.data[key]['AFBLM'] = AFBLMdict[key]['full'].copy()
    
    # Save object with Pickle
    with open(f'OCS_{key}_AFBLM_220722.pickle', 'wb') as handle:
        print(f'Saving {key} with pickle')
        pickle.dump(AFBLMdict[key]['full'], handle, protocol=pickle.HIGHEST_PROTOCOL)
    
# Save full class with Pickle
# NOW GIVING ERRORS ON RERUN (although worked previously), AttributeError: Can't pickle local object 'plotTypeSelector.<locals>.<lambda>'
# NOT SURE WHERE THIS IS COMING FROM AS YET
# UPDATE: IN data.data['ADM']['plots']['ADM']['pData'].attrs['pTypeDetails']
# TODO: GENERAL DEBUG FOR ACCIDENTAL FUNCTION PASSING HERE!!!!
# import pickle
with open(f'OCS_AFBLM_220722.pickle', 'wb') as handle:
    data.data['ADM']['plots']['ADM']['pData'].attrs['pTypeDetails'] = []  # UGH, this breaks pickle unless removed.
    pickle.dump(data, handle, protocol=pickle.HIGHEST_PROTOCOL)
    print(f'Saved data (full class) with pickle')